% 这个模型演示了增透膜的基本原理
% 波会在膜的前后表面反射。
% 只要设置膜的合适厚度，让前后表面的反射波的波程差等于波长的一半，那么这两个反射波就会相互抵消，从而减小反射
% 2d = lambda/2, d = lambda/4，膜的厚度至少是波长的1/4.
% *假定波速：入射介质>膜>出射介质。如果波速的顺序不同，你可能需要重新计算膜的厚度。
% CC-BY 4

clc
clear

n=2000;
dx = 0.01;
dt = 0.001;
c=[2;1.5;1]; %波速，分别是入射介质，膜，出射介质
T=1;
w=2*pi/T;

u0 = zeros(n,1);
u1 = zeros(n,1);
u2 = zeros(n,1);

k=0;

imgindex = 0;

Areflection = (c(3)-c(1))/(c(1)+c(3)); % 橙色虚线表示没有膜时的反射、透射波振幅。
Atransmission = 2*c(3)/(c(1)+c(3));

figure()
for t = 0:dt:20
  for i = 2:n-1
    if i < n/2
      j=1;
    elseif i >= n/2 && i <= n/2 +c(2)*T/4/dx
      j=2;
    else
      j=3;
    endif
    u2(i) = 2*u1(i) - u0(i) + (c(j)*dt/dx).^2*(u1(i-1)-2*u1(i) + u1(i+1));
  endfor

  if t > 4*T
    u2(1) = 0;
  else
    u2(1) = sin(w*t);
  endif
  u2(n) = 0;

  u0 = u1;
  u1 = u2;

  if mod(k,50) == 1
    clf
    hold on
    axis equal
    axis([0 20 -2 2])

    _x = linspace(0,dx*(n-1),n);
    _x = _x';
    plot(_x,u2)

    line([n/2*dx, n/2*dx], [-2 2]) %interface
    line([n/2*dx+c(2)*T/4, n/2*dx+c(2)*T/4], [-2 2]) %interface

    line([0 n/2*dx],  [Areflection Areflection],'linestyle','--')
    line([0 n/2*dx],  [-Areflection -Areflection],'linestyle','--')
    line([n/2*dx n*dx],  [Atransmission Atransmission],'linestyle','--')
    line([n/2*dx n*dx],  [-Atransmission -Atransmission],'linestyle','--')
    drawnow;
    pause(0.01)

    %print(['wave_transmission' num2str(imgindex) '.png'],'-dpng');
    imgindex++;
  end
  k++;
end

